Comparing Least Squares Calculations

Many statistics methods require one or more least squares problems to be solved. There are several ways to perform this calculation, using objects from the base R system and using objects in the classes defined in the Matrix package. We compare the speed of some of these methods on a very small example and on a example for which the model matrix is large and sparse. 1 Linear least squares calculations Many statistical techniques require least squares solutions β̂ = argmin β ‖y −Xβ‖ 2 (1) where X is an n × p model matrix (p ≤ n), y is n-dimensional and β is p dimensional. Most statistics texts state that the solution to (1) is β̂ = ( XX ) −1 Xy (2) when X has full column rank (i.e. the columns of X are linearly independent) and all too frequently it is calculated in exactly this way. 1.1 A small example As an example, let’s create a model matrix, mm, and corresponding response vector, y, for a simple linear regression model using the Formaldehyde data. > data(Formaldehyde) > str(Formaldehyde) 'data.frame': 6 obs. of 2 variables: $ carb : num 0.1 0.3 0.5 0.6 0.7 0.9 $ optden: num 0.086 0.269 0.446 0.538 0.626 0.782